function printout_P_b_c_from_omega(CONSTS, plot_data)

    num_points = 100;
    omega_p = CONSTS.omega_p;
    omega_H = CONSTS.omega_H;

    f_vec = linspace(omega_H/(2*pi), omega_p/(2*pi), num_points)';
    P_c_vec = zeros(size(f_vec));
    P_b_vec = zeros(size(f_vec));

    for i = 1:size(f_vec, 1)
        f = f_vec(i);
        LOCAL_CONSTS = CONSTS;
        LOCAL_CONSTS.f = f;
        LOCAL_CONSTS = update_CONSTS(LOCAL_CONSTS);

        P_b = LOCAL_CONSTS.P_b;
        P_c = LOCAL_CONSTS.P_c;

        P_b_vec(i) = P_b;
        P_c_vec(i) = P_c;

    end
    
    if (plot_data)
        figure; plot((2*pi)*f_vec, real(P_b_vec), 'b.-', (2*pi)*f_vec, real(P_c_vec), 'r.-'); title('P_b &P_c');
                 xlabel('\omega'); ylabel('P_c, P_b'); legend('P_b','P_c', 'location', 'SouthEast');
    end

end